Potential energy curves of graphene on different metal surfaces from the random 

phase approximation 
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We apply the Random Phase Approximation (RPA) to calculate the potential energy curves 
for graphene adsorbed on the Cu(lll), Ni(lll), Co(0001), Pd(lll), Pt(lll), Ag(lll), Au(lll), 
and Al(lll) metal surfaces. The results are compared with potential energy curves calculated 
with different semi-local density approximations and a non-local van der Waals functional, ft is 
shown that only the RPA functional captures both the weak covalent interactions and dispersive 
forces which are equally important for these systems. The RPA energies are evaluated in a plane 
representation using the projector-augmented wave method and we document our implementation 
in the electronic structure code GPAW, by comparing results for molecules and solids with previous 
ab initio RPA calculations. 

PACS numbers: 71.15.Nc, 73.22.Pr, 81.05.ue 
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I. INTRODUCTION 

The RPA was first introduced by Bohm and Pinesi~— 
more than 60 years ago and thus predates both Kohn- 
Sham Density Functional Theory (DFT) and the for- 
mal developments in many-body perturbation theory^. 
A most important property of the RPA is the explicit 
incorporation of screening in correlated quantities. The 
screening allows one to treat the electrons in metals as 
nearly independent " quasi-particles" interacting through 
a screened effective Coulomb interaction and explains 
why the single-particle picture often gives a decent de- 
scription of solids, despite the large strength of the bare 
Coulomb interaction^ The RPA takes screening into ac- 
count by summing a certain class of Feynman diagrams to 
infinite order in the Coulomb interaction and this allows 
one to evaluate correlation energies of metallic systems, 
which diverge in perturbative treatments £ 

In the context of DFT, the correlation energy can 
be expressed in terms of the interacting response func- 
tion using the adiabatic-connection and fluctuation- 
dissipation theorem (ACFD)»1 Time-dependent DFT re- 
lates the interacting response function to the Kohn- 
Sham response function through an exchange-correlation 
kernel and RPA is then the simplest possible approxi- 
mation where the kernel is neglected all together^ In 
principle, it is possible to calculate a local RPA po- 
tential using the optimized effective potential method, 
and solve the resulting Kohn-Sham equations for a self- 
consistent RPA density and total energy. However, the 
huge computational cost of such calculations has so far 
limited the self-consistent approach to atoms&i^ and sim- 
ple molecules^ Instead, it is often assumed that the ef- 
fect of selfconsistency is of minor concern and RPA cal- 
culations are performed non-selfconsistently using Kohn- 
Sham or Hartree-Fock eigenstates and eigenenergies. The 
approach then becomes equivalent to RPA from pertur- 
bation theory. 

The use of RPA as a tool for ab initio total energy 



calculations was pioneered by Furch o 12 ' 13 who calcu- 
lated the RPA atomization energies of small molecules 
and found that RPA has a general tendency to under- 
bind. It was also demonstrated that the atomic limit 
of N 2 dissociation was reproduced by RPA if the ref- 
erence is taken with respect to the RPA energy of two 
isolated N atoms. Subsequently, the performance of 
RPA has been examined systematically for molecular 
dissociatio n 14 ' 15 , cohesive energies of solids^ - — surface 
properties and adsorbates^ barrier heights ; 14 ' 15 ioniza- 
tion potentials^ 3 - van der Waals bonded dimers,— and van 
der Waals bonded two-dimensional materials £2r22. Com- 
pared to DFT calculations, the computational load of 
RPA calculations can represent a significant barrier for 
applications to large electronic systems. Nevertheless, 
due to increasing access to high performance computa- 
tional resources, there is a rapidly growing interest in the 
method and RPA is now slowly emerging as a standard 
tool in the electronic structure community. 

In general, RPA seems to be inferior to Generalized 
Gradient Approximations (GGA) and hybrid exchange- 
correlation functionals for the description of covalent 
bonds. However, the non-local nature of RPA makes its 
superior to any semi-local or hybrid functional for dis- 
persive interactions. Several attempts have been made 
to construct effective non-local van der Waals function- 
als, which capture long range dispersive interactions and 
are comparable to GGA calculations in computational 
requirements i^ 3 - - — In many cases, such functionals have 
been successful, but the approach is rather sensitive to 
the choice of exchange kernel and typically fails to give a 
qualitative description if both covalent and dispersive in- 
teractions are important j 28 ' 29 In contrast, the RPA corre- 
lation energy is naturally combined with exact exchange, 
does not rely rely on error cancellation or any fitted pa- 
rameters, and are able to describe intricate bonds with 
mixed covalent and dispersive character . 30 ' 31 For an ac- 
curate description of both strong covalent bonds and dis- 
persive interactions, it is necessary to apply beyond-RPA 
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methods such as screened second order exchange^ time- 
dependent EXXj 33 or renormalized adiabatic kernels^ 
However, such developments are outside the scope of the 
present paper and we will focus on RPA in the following. 

Two-dimensional layered compounds such as graphite, 
hexagonal boron nitride, and transition metal dichalco- 
genides, constitute a particular class of materials where it 
is vital to incorporate dispersive interactions in order to 
obtain a quantitative description of the bulk properties. 
RPA has been shown to provide an accurate description 
of van der Waals bonds in these system s 20 ' 21 ' 35 , whereas 
semi-local and van der Waals functionals can give rise 
to qualitatively wrong results. Moreover, the discovery 
and characterization of isolated graphene sheets; 36 ' 37 has 
triggered a vast amount of research in this intriguing two- 
dimensional material. In particular, graphene shows a 
remarkably high intrinsic carrier mobility, and therefore 
seems very well suited for nanoscale electronics devices. 
For such applications, the coupling to metal contacts 
plays a fundamental role and measurements show that 
graphene binds very differently on various metal surfaces. 
Understanding the interactions between graphene and 
metal surfaces, therefore becomes a most important task 
since the adsorption geometry and bond distance may 
have drastic consequences for the electronic structure and 
transport properties of adsorbed graphene layers. For 
example, experiments have demonstrated that Pd(lll), 
Co(0001), and Ni(lll) induce a band gap in adsorbed 
graphene sheets, which signals a covalent bond with the 
metal 3 ^— . In contrast, adsorption on Cu(lll), Ag(lll), 
Au(lll), and Pt(lll) do not change the graphene band 
structure significantly 41 - - — . On the theoretical side, most 
studies have been limited to semi-local approximations 44 
and van der Waals functional s 45 ' 46 . While some agree- 
ment with experiment was obtained in Ref. 46] for a 
certain van der Waals functional, the large discrepancy 
between various choices of functional is clearly unsatis- 
factory. 

Here we apply RPA to calculate the binding en- 
ergy curves of graphene on Cu(lll), Ni(lll), Co(0001), 
Pd(lll), Pt(lll), Ag(lll), Au(lll), and Al(lll) metal 
surfaces. The results for Cu(lll), Ni(lll), and Co(0001) 
have been obtained previousl y 30 ' 31 but are reproduced 
here for completeness. We also show that the slight dis- 
crepancy between the RPA curve for graphene on Ni(lll) 
in Refs. (30j and l3l| were caused by insufficient /c-point 
sampling in Ref. [30j . For all the metal surfaces except 
Pd(lll), we find good agreement with experiment. The 
deviation in the case of Pd, is most likely related to the 
large discrepancy between the metal and graphene unit 
cells and a proper Moire structure is needed in order to 
compare with experiments in this case. 

The paper is organized as follows. In section |H] we out- 
line the general method used to obtain RPA total energies 
and present details on the plane wave implementation ap- 
plied in the present work. In section Hill the RPA poten- 
tial energy curves for graphene adsorbed on 8 different 
metal surfaces are presented and compared with semi- 



local approximations for the exchange-correlation energy 
and a standard van der Waals functional. We then asses 
the quality and wide applicability of the method and im- 
plementation by benchmarking calculated results for dis- 
sociation of graphite, properties of bulk solid state sys- 
tems and molecular atomization energies. In appendices 
I A II and I A 21 we present detailed convergence tests for the 
RPA potential energy curves of graphene on Ni(lll) and 
for the atomization energy of the CO molecule. 



II. METHOD 
A. Theory 

Using the adiabatic connection and fluctuation- 
dissipation theorem (ACFD), the exchange-correlation 
energy can be written as: 

E xc = - J d\J°° ^-Tr{v[n2nS(u}) + X A M]}, (1) 

where n(r, r') = n(v)5(r — r') and v is the Coulomb in- 
teraction. Here n(r) is the density, which by definition 
is constant along the adiabatic connection and is 
the interacting response function of a system with v — > Xv 
evaluated at imaginary frequencies. It is standard prac- 
tice to divide E xc into an exchange part E x obtained by 
setting A = in the integrand and a correlation part E c , 
which is the remainder. One then obtains 

E c = -J dX J^^Tr{v[ X X (iu>)-x KS (w)}}, (3) 

where is the response function of the non- 

interacting Kohn-Sham system. A major advantage of 
this separation is that the exchange energy can be eval- 
uated exactly and one only needs to approximate % A to 
obtain E c . 

The Random phase approximation for the interacting 
response function can be derived in several ways, but in 
the present context it is convenient to use time-dependent 
density functional theory, from which it is it is straight- 
forward to show that 

(4) 

where f£ c is the exchange-correlation kernel. The RPA 
is then obtained by taking f£ c = and one is left with 

= [! - X KS (^)Xv] -\ KS (iu>). (5) 

Inserting this into the expression for the correlation en- 
ergy and carrying out the coupling constant integration 
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yields 

E RPA 



dw r 
2^ 



-Tr 



{\n[l~v X KS (^)}+v X KS ^)y 



(6) 



For spin polarized systems the correlation energy involves 
the spin summed response function x x — ^2 aa i Xaa 1 - Us- 
ing that v is independent of spin, it is straightforward 
to show that x A satisfies Eq. ([5]) if \ KS replaced by 
xf S + X^ S ■ This would not be true if a spin-dependent 
f xc were included in Eq. (T5|) and comprises a major sim- 
plification of RPA calculations involving spin polarized 
systems. 



B. Plane wave implementation 

For solid state systems it is convenient to evaluate the 
response function in a plane wave representation. The 
number of plane waves required at a given energy cutoff 
scales as Nq ~ V ce u, which determines the dimension of 
the response function. For isolated atoms and molecules 
where large unit cells is needed in order to decouple pe- 
riodic images, the response function may become pro- 
hibitly large and the implementation is not well suited 
for large molecular systems. However, as will be shown 
below, it is possible to calculate atomization energies for 
small molecules although the computational time is much 
larger than implementations utilizing atomic basis sets. 

In a plane wave basis the Kohn-Sham response func- 
tion is 

KS I ■ \ 1 fnk ~ /n'k+q 

Xgg'(<L««0 = Tr 1^ 1^ ~T 7 

x(V„k|e- < ( cl + G )- r |^k +q }<V'n'k +q |e i ( q+G '> r |Vnk> ) 

where e„k are Kohn-Sham eigenvalues, / n k are occupa- 
tion numbers, and \ipnk) are the Kohn-Sham eigenstates 
normalized in the unit cell with volume V. The trace 
in Eq. © then becomes a trace over plane waves and 
a Brillouin Zone integral over q which is sampled on a 
uniform grid: 



E. 



RPA 



duj 1 
2^W a 



■E 

qEBZ 



(8) 



x Tr 



{ ln[l - «(q)x KS (q, iu)] + v(q) X KS (<l, «•>)}• 

The plane wave representation of the coulomb interac- 
tion is v := 47T(5gg'/I < 1 + G| 2 and the trace of the loga- 
rithm is most easily evaluated by using that Tr[ln(yl)] = 
ln[dct(yl)]. The exact exchange energy Eq. ([2]) becomes 



E 



EXX 



iVqiVk 



n,n' k,qG.B2 



x ^ VG (q)|(^k|e-^ +G > r |^k +q )| 2 (9) 



The expression is derived from the ACFD and differs 
from the standard expression for exact exchange energy, 
if the occupation numbers are not integer valued. How- 
ever, as discussed in Ref. [TJ], it is natural to apply Eq. 
(|9]) when the exact exchange energy is combined with 
the RPA correlation energy. For metals, it is customary 
to aid convergence by smearing the occupation factors by 
an artificial electronic temperature and it has been shown 
empirically that Eq. (|9|) is less sensitive to the width of 
the artificial smearing function / n k than the standard 
expression for exact exchange. 

The q = terms in Eqs. © and © require a special 
treatment since t?(q) diverges as q — > when G = 0. The 
divergence is, however, integrable and the terms yield a 
finite contribution. For the exact exchange part we ap- 
ply the method of Gygy and Balderesch i 47 ' 48 where the 
Coulomb interaction is multiplied by a Gaussian regular- 
ization and the q = term can be integrated analytically 
in the limit of infinitely dense Appoint sampling. The 
correlation energy is evaluated by replacing |i/Vik+q) by 
its first order perturbative expansion in the k • p Bloch 
Hamiltonian. For n ^ n' this yields 

(V>nk|e 4 |VVk+q)q^0 = (10) 

^n'k ^nk 

and the q on the right hand side cancels the diverg- 
ing Coulomb interaction in Eq. ([8]). The limit clearly 
depends on the polarization of q and for non-isotropic 
systems we evaluate the q = contribution to the cor- 
relation energy by averaging over non-equivalent polar- 
izations. The method may fail for systems with a high 
density of states near the Fermi level such as certain tran- 
sition metals, since the denominator in Eq. (|10[) approach 
zero for low energy transitions. In principle, the problem 
should be solved by using degenerate perturbation the- 
ory, but there is no unique way of defining which states 
to treat as degenerate for a given k. Instead we take a 
pragmatic point of view and simply exclude the q = 
term in the evaluation of (|5]) and © for systems involv- 
ing transition metals. This procedure has been shown 
to exhibit fast convergence for E x + E c with respect to 
fc-point sampling 1 ^, whereas the individual exchange and 
correlation terms converge rather slowly when q = is 
excluded. 

It should be noted that the g-point grid is completely 
determined by the fc-point sampling of the initial self- 
consistent Kohn-Sham calculation, since q has to cor- 
respond to a fc-point difference. For all calculations we 
use Gamma-centered uniform fc-point grids since the fc- 
points and g-points then coincide. This ensures a much 
more efficient symmetry reduction of the g-points than if 
a shifted fc-point grid were to be used. To evaluate the re- 
sponse function we usually choose a cutoff energy E* ut , 
which is smaller than the cutoff used to obtain the in- 
put eigenstates and eigenenergies and the set of included 
plane waves are determined by |q+G| 2 /2 < E* ut . The di- 
mension of the response function xgg' (q) thus depends 
on q, which ensures a smooth dependence on the cut- 



4 



off energy for periodic systems. In all calculations, the 
number of bands used to evaluate the response function 
is set equal to the number of plane waves determined by 
the cutoff. This approach is appealing, since the RPA 
calculations then only depends on a single convergence 
parameter. Furthermore, it will be more straightforward 
to compare with properties of the homogeneous electron 
gas and the Lindhard function, since in that case, the 
eigenstates coincide with plane waves. 

The response function is evaluated at the imaginary 
frequency axis, where it varies rather smoothly and al- 
lows for an efficient numerical integration. Typically, the 
density of states near the Fermi level determines how 
much structure the response function exhibits near u — 
0. The frequency integration in Eq. (JHJ is carried out us- 
ing 16 Gauss-Legendre points with a weight function en- 
suring that the integral of f(x) oc a^ 1 /- 8 " 1 ) exp [—ax 1 / 5 ] 
is reproduced exactly Here a is determined by the high- 
est frequency point, which we position at 800 eV for 
all calculations. B determines the density of frequency 
points close to ui — and we use B = 2.0 for systems 
with a gap and B = 2.5 for metals. With this frequency 
sampling the RPA correlation energies are converged to 
within a few meV. 

Since, the present approach is not self-consistent, one 
has to choose a set of orbitals, on which Eqs. ©-© 
are evaluated. We have compared the RPA potential en- 
ergy surface for graphene on Ni(lll) using self-consistent 
PBE orbitals with the result obtained with self-consistent 
LDA orbitals and the results are very similar. In other 
cases, however, there may be a significant dependence 
on the initial orbitals. For example, the RPA atomiza- 
tion energy of O2 and CO was shown to differ by ~ 0.4 
eV when comparing LDA and PBE initial orbitals 34 . 
Even larger differences have been observed when com- 
paring Hartree-Fock and PBE initial orbitals 4 — and the 
most accurate RPA results are obtained when combining 
self-consistent EXX with RPA correlation energies eval- 
uated on PBE or similar initial orbitals. Alternatively, 
the effect of non-selfconsistency can be corrected by in- 
cluding single-excitation contributions to the correlation 
energy 49 . Unless stated otherwise, all RPA calculations 
below are performed with PBE Kohn-Sham orbitals and 
eigenvalues. 

The response function, EXX, and RPA expressions 
Eqs. d)-© has been implemented in GPAW— ~— , which 
is a Density Functional Theory code based on the pro- 
jector augmented wave (PAW) method^ 3 -. We refer to 
Ref. [54| for details on the PAW implementation of the 
response function. 



III. RESULTS 

A. Graphene on metal surfaces 

The main result of the present paper is the RPA calcu- 
lation of potential energy surfaces of graphene on 8 dif- 





FIG. 1: Left: Minimal unit cell used for Co(0001), Cu(lll) 
and Ni(lll). Right: y/3 x unit cell used for Pd(lll), 
Pt(lll), Au(lll), Ag(lll), and Al(lll). 



Metal: 


Ni 


Co 


Cu 


Pd 


Pt 


Au 


Ag 


Al 


a (A) 
AE (meV) 


2.49 
6 


2.51 
21 


2.56 
92 


2.38 
91 


2.40 
52 


2.50 
13 


2.51 
21 


2.48 
2 



TABLE I: Value of the graphene lattice parameter a when 
matched to the various metal surfaces with experimental lat- 
tice parameters. We also display the energy per C atom AE 
(calculated with the PBE functional) needed to stretch an iso- 
lated graphene sheet to the metallic lattice parameter. The 
experimental lattice parameter of isolated graphene is 2.46 A. 



ferent metal surfaces. The calculation has already been 
carried out for N^lll) 3 ^, Cu(lll) and 00(0001)22 
which all can be done with the minimal lxl surface 
unit cell. Here we extend the calculation to include 
Pd(lll), Pt(lll), Au(lll), Ag(lll), and Al(lll), where 
the y/3 x v3 surface unit cell is needed in order to obtain 
a periodic system which is compatible with the graphene 
lattice distance. In all calculations we have used the 
experimental metal lattice parameter, and stretched or 
squeezed graphene to match the unit cell. In Table Q] we 
present applied lattice parameters along with the energy 
required to stretch an isolated graphene sheet to match 
the lattice. For Ni, Co, and Cu we have a = d and for 
Ag, Au, Pd, Pt, and Al we have a — 2d/y/3 for a metal 
nearest neighbor distance d. In the case of Cu, Pd and 
Pt the deviation from the applied unit cell is particularly 
bad and it is expected that a more complicated Moire 
pattern is needed in order for the graphene and metal 
surface to come in registry. Nevertheless, it is interesting 
to compare the performance of different functionals even 
though these structures are not observed experimentally. 

For all calculations except the RPA correlation energy, 
we used a plane wave cutoff of 600 eV. For the RPA 
correlation energy we used the orbitals and eigenvalues 
obtained with 600 eV cutoff and evaluated the response 
function at a cutoff of 200 eV for the small systems and 
150 eV for the large systems. The number of bands in- 
cluded in the response function were set equal to the 
number of plane waves defined by the cutoff energy. The 
metal surface was simulated using four atomic layers and 
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the repeated images were separated by 20 A of vacuum. 
For the Ni(lll) and Co(0001) slabs the calculations were 
spin-polarized. A 16 x 16 gamma-centered k-point mesh 
was used for Ni and Co whereas an 12 x 12 grid was used 
for Cu and 6x6 grids were used for the large systems. 
We return to the issue of k-point and cutoff convergence 
in appendix I A 11 

In Fig. [2] we show the potential energy curves of 
graphene on the Ni(lll), Co(0001), Cu(lll), Au(lll), 
Ag(lll), Pt(lll), Pd(lll), and Al(lll) calculated with 
LDA, PBE, RPBE, a van der Waals functional^ (vdW- 
DF) , and RPA. The binding energies and equilibrium dis- 
tances are summarizes in Tabs. ITTIand llTll LDA predicts 
strong binding (~ 70 — 260 meV) and small binding dis- 
tance (2.0-2.21 A) for Ni(lll), Co(0001) and Cu(lll) 
and weak binding at ~ 3.3 A for the rest of the metals. 
It should be noted that if the LDA optimized lattice con- 
stant is used for Cu the binding is weaker and similar to 
the LDA curve of Pd(lll)4^ PBE predicts a very weak 
bond at ~ 4.4 A except in the case of Co(0001) where a 
minimum close to the surface (~ 2.0 A) appears. This 
feature is also observed for Ni(lll) where a local min- 
imum appears close to the surface, but in this case it 
is unstable with respect to the desorbed graphene. The 
RPBE functional predicts very weak binding far from the 
surface (~ 4.4 A) for all the systems. The vdW-DF also 
gives very similar results for all the metals with an equi- 
librium distance of ~ 3.75 A and a binding energy of 
- 40 meV. For Ni(lll) and Co(0001) RPA produces two 
distinctive minima at ~ 2.2 A and ~ 3.25 A and in both 
cases the global minimum is the one close to the surface. 
For the rest of the systems RPA predicts an equilibrium 
distance at ~ 3.3 A, but with much larger binding ener- 
gies than any of the other functionals. 

It has previously been shown that the electronic struc- 
ture of graphene adsorbed on Ni(lll) and Co(0001) 
is severely modified at the minimum close to the sur- 
face, whereas the graphene electronic states do not hy- 
bridize with the metallic states at ~ 3.25 A*22, Both 
the small binding distances and the modified electronic 
structure at these distances are in good agreement with 



experiments 



38,40 



For the remaining systems a direct 



comparison with experiments is not possible, since ex- 
tended Moire patterns are observed due to a mismatch 
of lattice parameters. Nevertheless, we can use these 
systems as a test set for comparing the performance of 
different functionals. RPA and vdW-DF is the only non- 
local functionals considered and they both capture the 
slowly decaying long distance tail originating from dis- 
persive interactions. However, at intermediate distances 
where both dispersive and covalent interactions are im- 
portant the two functionals deviate significantly. The 
vdW-DF gives very similar results for all systems and 
does not seem to capture the differences in surface elec- 
tronic structure when the surface is approached. In fact, 
changing the local part of the vdW-DF has been shown 
to give rise to qualitatively different energy curvesi 29 ! 31 
It thus seems difficult to choose an accurate vdW-DF for 



Metal: 


Ni 


Co 


Cu 


Pd 


Pt 


Au 


Ag 


Al 


E^ DA (meV) 


188 


259 


72 


43 


36 


34 


30 


29 


E PBE (meV) 


2 


29 


2 


4 


5 


2 


2 


2 


EW BE (meV) 


1 


1 


1 


2 


3 


1 


1 


1 


E v B dw (meV) 


39 


38 


39 


40 


42 


40 


36 


36 


Eg PA (meV) 


70 


78 


68 


90 


84 


95 


78 


52 



TABLE II: Binding energies per C atom at the equilibrium 
distance to the surface calculated with different functionals. 



Metal: 


Ni 


Co 


Cu 


Pd 


Pt 


Au 


Ag 


Al 


d L B DA (A) 


2.00 


2.01 


2.21 


3.00 


3.35 


3.32 


3.22 


3.44 


d p B BE (A) 


4.33 


2.12 


4.33 


4.25 


4.40 


4.53 


4.47 


4.55 


d% PBE (A) 


5.48 


5.60 


5.54 


5.26 


5.22 


5.43 


5.57 


5.61 


d v B dw (A) 


3.73 


3.80 


3.80 


3.73 


3.84 


3.82 


3.84 


3.99 


d R B PA (A) 


2.19 


2.27 


3.09 


3.34 


3.42 


3.22 


3.31 


3.51 



TABLE III: Equilibrium distances to the surface calculated 
with different functionals. 



a given system of this type, whereas RPA constitutes a 
unique functional that do not involve arbitrary choices 
for exchange and (local) correlation. 

We note that one would expect the RPA energy curves 
at long distances to be well described by vdW-DF func- 
tionals since the local contributions to exchange and cor- 
relation then vanishes. This seems to be case for all 
the metals except Pd(lll) and Pt(lll). It is interest- 
ing that for these two metals the pure HF energy curves 
produce weak minima at 5 A with binding energies of 5 
and 8 meV respectively. For all other surfaces the HF 
energy curves are purely repulsive in this region. The ex- 
change functional used in the vdW-DF considered here 
is the revPBE. 55 This is very similar to the RPBE func- 
tional, which gives completely similar structure at large 
distances from the surface. It is thus very likely that the 
deviations between RPA and vdW-DF at a distance of 
~ 5 — 6 A from the surface is due to small local exchange- 
correlation effects, which are not well described by the 
present vdW-DF. The fact that (semi) local exchange- 
correlation effects are important at distances of ~ 5 A for 
the surface is also supported by recent calculations with 
the M06-L functional, 5e - which accurately reproduce the 
RPA energy curve for Ni(lll)££ 



B. Graphite 

A very important accomplishment of the RPA method 
is the demonstration of a correct description of the cohe- 
sive properties of graphite^ 1 RPA gives excellent agree- 
ment with the experimental interlayer distance and inter- 
layer binding energy which so far only has been calculated 
accurately with quantum Monte Carlo methods^ Here 
we reproduce the main results of Ref. 21( in order to 
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asses the performance of the present implementation. 

For the graphite calculations we used a Gamma cen- 
tered 26 x 26 x 8 fc-point grid for the DFT and HF calcu- 
lations and a 14 x 14 x 6 fc-point for the RPA calculations. 
We have used a plane wave energy cutoff of 800 eV for 



the HF and DFT calculations. To obtain the RPA in- 
terlayer binding energy of graphite, one has to compare 
correlation energies evaluated in different unit cells corre- 
sponding to different values of the interlayer separation d. 
Therefore, one cannot use a single low energy cutoff and 
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-5.8 

> 








E cut [eV] 

FIG. 3: RPA correlation energy per atom as a function of 
cutoff energy at d = 3.34 A. For cutoff energies above f75 
eV the extrapolation scheme is rather accurate and the un- 
certainty on the extrapolated result is on the order of 5 meV 
and f meV for two-point extrapolation and linear regression 
respectively. 




3.0 3.5 4.0 4.5 5.0 5.5 6.0 
d [A] 



FIG. 4: Potential energy surfaces for graphite obtained with 
RPA, LDA, PBE and a van der Waals functional. The error 
bars on the RPA calculations are obtained from the linear 
regression applied to calculate the extrapolated correlation 
energies (Eq. ((TTJ) 



rely on error cancellation between energy differences as 
e.g. in the case of graphene on metal surfaces. Instead, 
for a given value of d one has to obtain the converged 
value of the correlation energy corresponding to infinite 
cutoff energy. As shown empirically in Ref. (l6| . the 
correlation energy at high values of the cutoff scales as 

3/2 

~ E cut and the converged RPA correlation energy can 
be obtained by fitting the function 



E^ A {E cut ) = E^ A + 



.4 



E, 



3/2 • 



(11) 



to a sequence of cutoff energies and the associated cor- 
relation energies. In Fig. [3] we show the extrapola- 
tion at the equilibrium distance of d = 3.34 A. The fit 
to Eq. (fTTj) appears rather accurate and the extrapo- 
lated correlation energy ranges from —6.958 eV to —6.949 
eV when any two points between E cut — 175 eV and 
Ecut — 275 eV are used. However, to obtain a mean- 
ingful binding energy curve, an accuracy of ~ 1 meV is 
needed and a two-point extrapolation is not sufficient. 
Instead, we perform linear regression on the five points: 
E cut e {175,200,225,250,275} resulting in the extrap- 
olated correlation energy E RPA = 6.955 ± 0.00095 eV. 

This extrapolation procedure is repeated for a range 
of different values of d and the result is shown in Fig. 0] 
along with the results obtained with the LDA, PBE, and 
a van der Waals functional^ The result is well known: 
LDA predicts a fortuitous equilibrium distance, which is 
in good agreement with experiment and the PBE curve is 
purely repulsive. The van der Waals functional seems to 
capture part of the dispersive interactions and predicts 
a larger binding energy than the semi-local functionals, 
however, at the wrong equilibrium distance. The RPA 
method gives a binding energy, which is in very good 



agreement with experiments and quantum Monte Carlo 
simulations^ and the correct equilibrium distance at d = 
3.34 A. Note that we obtain a slightly larger RPA binding 
energy per C atom (62 meV) than Lebegue et al<21 (48 
meV). The reason for this could be related to the use of 
different PAW setups for C. It is also possible that our 
/c-point sampling is not completely converged, since if we 
use the same fc-point sampling (14 x 14 x 6) for both 
Hartree-Fock and RPA correlation we get a binding of 
47 meV per Carbon atom. However, with this fc-point 
sampling the large distance tail of the Hartree-Fock PES 
is not converged and the total PES exhibits a spurious 
maximum at d ~ 5 A. 



C. Cohesive energies and lattice constants of solids 



In Refs. [HT3 the RPA method was demonstrated to 
yield bulk lattice constants in very good agreement with 
experiment and cohesive energies somewhat worse than 
the PBE functional. The absolute RPA correlation en- 
ergy is typically overestimated by 25-50 % 13 ' 15 for atoms 
and molecules and slightly less for solids, but RPA en- 
ergy differences are reproduced accurately when systems 
with similar electronic structure is compared. Lattice 
constants are determined by comparing very similar elec- 
tronic systems and are therefore well reproduced by RPA 
calculations. In constrast, the computation of cohesive 
energies requires comparison of atoms in the solid phase 
with the isolated atoms, which have a completely differ- 
ent electronic structure and RPA performs poorly in this 
case. 

In Tab. |lV]we display the cohesive energies of a selec- 
tion of solids calculated with PBE, EXX, and RPA. In 
all calculations we used a gamma-centered fc-point sam- 
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PBE 


EXX 


RPA 


Expt. 




@>exp &opt 


Q>exp &opt 


&exp @>opt 




c 


7.73 (7.72) 


5.16 (5.18) 


6.99 (7.00) 


7.55 


Si 


4.55 (4.55) 


2.82 (2.82) 


4.32 (4.39) 


4.68 


SiC 


6.38 (6.40) 


4.32 (4.36) 


5.96 (6.04) 


6.48 


Ge 


3.72 (3.71) 


2.05 (1.95) 


3.56 (3.59) 


3.92 


MgO 


4.97 (4.98) 


3.35 (3.47) 


4.85 (4.91) 


5.20 


Na 


1.08 (1.08) 


0.20 (0.23) 


0.98 (1.00) 


1.12 


Pd 


3.68 (3.74) 


-1.44 (-1.26) 


3.51 (3.41) 


3.94 


Rh 


5.61 (5.74) 


-3.01 (-2.88) 


5.10 (5.05) 


5.78 


Cu 


3.40 (3.48) 


-0.23 (0.03) 


3.20 (3.36) 


3.52 


MAE 


0.16 (0.13) 


3.22 (3.14) 


0.41 (0.38) 







PBE 


EXX 


RPA 


Expt. 


C 


3.57 (3.57) 


3.55 (3.54) 


3.57 (3.57) 


3.55 


Si 


5.48 (5.47) 


5.49 (5.48) 


5.45 (5.43) 


5.42 


Na 


4.20 (4.20) 


4.47 (4.49) 


4.29 (4.18) 


4.21 


Pd 


3.95 (3.94) 


4.03 (4.00) 


3.90 (3.90) 


3.88 



TABLE V: Optimal lattice constants of a few solids. Experi- 
mental values are corrected for zero-point anharmonic effects. 
All numbers are in A. Numbers in brackets are taken from 
Ref. [3. 



D. Dissociation of molecules 



TABLE IV: Cohesive energies of solids evaluated at the exper- 
imental lattice constant corrected for zero-point anharmonic 
effects. Numbers in brackets are at optimized lattice constant 
and are taken from Ref. pl| • Experimental cohesive energies 
are corrected for zero point energy. All numbers are in eV. 



pling of 12 x 12 x 12 for the solid. For the calculation of 
the isolated atoms, the periodic images were separated 
by 8 A except Na, where a separation of 10 A was used. 
The RPA energy differences were calculated at different 
cutoff energies and extrapolated to infinity. A two-point 
extrapolation using Eq. (fTT|) with either {250,300} eV 
or {350,400} eV, yielded results differing by ~ 2 meV 
(Si, Ge, Na) to 50 meV (Pd, CuL We find good agree- 
ment with the results of Ref. [18| with a deviation of 
0.01 — 0.1 eV. It should be noted that the results of 
Ref. [l8j were calculated at optimized lattice constants 
whereas the present results are at the experimental lattice 
constants. This is part of the reason why our calculated 
cohesive energies are generally smaller than those of Ref. 
18] and the largest deviation is seen for EXX applied 
to metals where a large difference from experimental lat- 
tice constants is observed. We should also remark that 
the PAW setups used in the present work have not been 
optimized for RPA calculations as in Ref. [l~8T |. 

In Tab. [V] we show the calculated lattice constants 
of C, Si, Na, and Pd. Again we find good agreement 
with the results of Ref. [18]. The results were obtained 
by calculating 7 — 9 energy points in the vicinity of the 
experimental lattice constant and fitting a third order 
inverse polynomial to the energy- volume curve i^2r— The 
RPA results were obtained by a two-point extrapolation 
with {250,300} eV using Eq. ((TTJ- In principle, this 
approach also gives the bulk modulus as the curvature 
at the minimum. However, for the RPA calculations, 
the present extrapolation scheme is not accurate enough 
for this purpose. Linear regression involving more cutoff 
points would be needed in order to produce a reliably 
RPA bulk modulus. 



The calculation of molecular atomization energies is 
not well suited for a plane wave implementation, since 
the number of plane waves included at a given energy 
cutoff, scales as the cube of the super cell size. There- 
fore, the dimension of the response function xgc quickly 
becomes prohibitly large when the super cell is increased 
and it becomes very difficult to compute RPA correlation 
energies is a plane wave basis. Nevertheless, we can cal- 
culate RPA atomization energies for small molecules and 
compare our implementation with codes using atomic ba- 
sis sets. 

To obtain the RPA correlation part of atomization 
energies, we use Eq. (fTTj) to extrapolate calculations 
performed at E cut <E {150,200,250,300,350,400} eV. 
While the absolute RPA correlation energies are hard to 
converge, extrapolated energy differences are converged 
when points in the range E cut € {300, 350, 400} are used. 
Thus, if the same unit cell is used for the calculation of 
the molecular correlation energy and the atomic correla- 
tion energies, the energy difference can be obtained by 
a two-point extrapolation using either E cut € {350,400} 
eV or E cut e {300, 350} eV. The extrapolated energy 
differences differ by at most 20 meV depending on which 
two cutoff energies are used. In contrast, the individ- 
ual correlation energies of atoms and molecules are much 
harder to converge and linear regression is needed in or- 
der to get a reliable extrapolated result. In general, larger 
unit cells tend to improve the accuracy of the extrapola- 
tion, since the larger number of plane waves results in a 
smoother cutoff dependence. In Appendix IA 21 we show 
various convergence test for the correlation energy of the 
CO molecule. 



1. Atomization energies 

We have computed the atomization energies of 12 small 
molecules and compared with the results obtained by 
Furch o 2 ' 62 using an atomic basis set approach. The cal- 
culations were performed on experimental geometries and 
the experimental atomization energy has been corrected 
for zero-point vibrational energies. For P and CI the cal- 
culations were performed in a supercell where the nearest 
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LDA 


PBE 


RPA - This work 


RPA - Ref. [62] 


Exp. 


H 2 


113 


105 


109 


109 


109 


N 2 


268 


244 


224 


223 


229 


2 


174 


144 


112 


113 


121 


CO 


299 


269 


243 


244 


259 


F 2 


78 


53 


30 


31 


39 


HF 


161 


142 


131 


133 


141 


H2O 


266 


234 


222 


224 


232 


C2H2 


460 


415 


383 


381 


405 


CH4 


462 


420 


404 


405 


419 


NH 3 


337 


302 


290 


290 


297 


Cl 2 


81 a 


65 


49 


50 


58 


P 2 


143 


121 


116 


116 


117 


MAE 


36 


8 


9 


9 







PBE 


PBEx 


HFQPBE 


RPA@PBE 


Ref. [64] 




2.86 


1.86 


-1.85 


2.48 


2.68 


a 3 n 


2.63 


2.06 


0.92 


2.25 


2.48 



a We were not able to converge the LDA energy of the isolated CI 
atom and the non-selfconsistent LDA energy evaluated at the PBE 
density was used here. 

TABLE VI: Atomization energies of small molecules. The 
results from Ref. [12] were performed with atomic orbital 
basis set. All RPA energies were performed with selfconsistent 
PBE orbitals and eigenenergies. All numbers are in kcal/mol 
(1 kcal/mol=43 meV). 



neighbor atoms of periodic images were separated by 8 A. 
For the rest of the elements a separation of 6 A was suf- 
ficient. The results are shown in Tab. I VII and we observe 
a close agreement with the results of Furche. 

The comparison with experimental values and the PBE 
functional is well known. RPA systematically underes- 
timates atomization energies and performs slight worse 
than the PBE functional, but significantly better than 
LDA. 



2. Static correlation of MgO dimer 

An accurate description of the MgO dimer ground 
state energy, represents a challenge for any single refer- 
ence ab initio method due to the multi reference nature 
of the ground state. 63 i 64 The two lowest lying electronic 
states are the singlet X 1 ^" 1 " and triplet a 3 n with the for- 
mer being favored by 0.2 eV. Both of these states have an 
"open shell" ionic character with Mg donating an elec- 
tron to O and are correlated with ionic diabatic dissoci- 
ation limits^ 3 - This results in a significant hybridization 
between valence and Rydberg states and the dimer in 
its ground states is not well characterized by a single 
Slater determinant. Moreover, while the adiabatic po- 
tential energy curve for the triplet naturally dissociates 
into the lowest lying atomic configuration 1 Mg+ 3 O, the 
singlet will dissociate into 1 Mg + 1 O. Here we exam- 
ine how the the atomization energy of the lowest sin- 
glet and triplet states are described by Hartree-Fock 
and RPA. Thus we calculate the atomization energies 



TABLE VII: Atomization energies of the lowest singlet and 
triplet states of the MgO dimer. All numbers are in eV. 



E a = Ef/s — EMg — Eo, where E t / S are the ground state 
energies of the triplet /singlet, 2?Mg is the energy of a sin- 
gle Mg atom in its singlet state, and Eo is the energy of 
a single O atom in the triplet state. 

In Table IVlTI we show the calculated atomization ener- 
gies of MgO using the PBE functional, Hartree-Fock, and 
RPA and compare with high-level correlated methods. 64 
The calculations were performed with fixed equilibrium 
geometries taken from Ref. 64] . PBE overestimates the 
atomization energies slightly, but predicts the correct or- 
der of adiabatic states. In contrast, HF is not able to 
capture the static correlation originating from the ionic 
configuration of the dimer and predicts the singlet to be 
unstable. Remarkably, RPA produces correlation ener- 
gies for the singlet and triplet, which differ by ~ 3.0 eV, 
but they correct the HF energies just right, such that 
the order of adiabatic states is restored. The total RPA 
atomization energies slightly underestimate the exact at- 
omization energies, which is in line with the trend previ- 
ously observed for small molecules. 

The reason for very different contributions from HF 
and RPA correlation, despite similar total energies, is 
most likely related to the second terms of Eqs. ([5} and 
From the point of view of the adiabatic connection, 
the separation into exchange and correlation is somewhat 
arbitrary, since x KS nas been added and subtracted from 
Eq. ([1]) . Clearly, if \ KS gives a large contribution to Eq. 
([2]), one would not expect the HF energy to be accurate. 
In contrast, when the RPA correlation energy is added, 
the contribution from \ KS cancels out and one is left 
with the two terms in Eq. ([1}. In this respect, the most 
surprising result in table IVHI is perhaps the fact that the 
PBE energies are so close to the exact result. This is a 
manifestation of the accurate error cancellation between 
exchange and correlation in the PBE functional. The 
PBE exchange energy is very far from the exact exchange 
(here represented by non-selfconsistent HF), but when 
correlation is included the PBE and RPA approach yield 
similar results with the correct ordering of states. 



3. The atomic limit of molecular dissociation 

A surprising feature of the RPA is the cor- 
rect description of the atomic limit of molecular 
dissociation i 12 ' 13 i 15 ' 65 Apparently, the non-perturbative 
nature of RPA captures the strong static correlation aris- 
ing in the atomic limit, which is a remarkable property of 
a single reference method. For example, semi-local DFT, 
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Hartree-Fock and coupled cluster typically yield dissocia- 
tion limits, which have to high energies. ~ However, RPA 
fails dramatically in the case of dissociation, which 
is completely free of electronic correlation. In fact, the 
RPA total energy of a single H atom is ~ —0.6 eV and 
the atomic limit of H2 dissociation only comes out right 
if taken with respect to the RPA reference of a single H 
atom. 

Here we will attempt to reproduce the well known en- 
ergy curves of H2 dissociation using our plane wave im- 
plementation. Of course, the atomic limit of molecu- 
lar dissociation is extremely difficult to reproduce with 
plane waves and periodic boundary conditions due to the 
large unit cells required. However, it is a nice test of the 
present implementation to see if the static correlation 
can be captured using plane waves and periodic bound- 
ary conditions. In the case of H2, convergence is very 
fast with respect to cutoff and we can manage to ob- 
tain converged RPA correlation energies using extrapo- 
lation with E cut e {100,150} eV. In Fig. [5] we show 
the dissociation curve of H 2 where we used a unit cell 
size of 2.5d x 2.5d x 3.5d with d being the H-H distance. 
We were not able to perform calculations with unit cell 
sizes beyond 12 x 12 x 16 A, but the trend of the RPA 
curve seems to agree very well with the results of Refs. 

However, note the results are not completely 



converged with respect to unit cell size. The spurious 
maximum at d = 3.7 A, also found in previous studies, is 
situated at E max = 0.66 eV. If we instead use unit cells 
of 2d x 2d x 3d and 3d x 3d x 3d we obtain E max = 0.56 
eV and E max = 0.72 eV respectively. It should be noted 
that the energy in Fig. [5] is taken with respect to two 
isolated H atoms for which RPA gives E^ PA = —0.57 
eV. On an absolute scale RPA would thus underestimate 
the entire energy curve by ~ 1.1 eV. 

In Fig. [5] we show the dissociation curve of H^ . Again, 
we emphasize that these molecular systems are far from 
the periodic systems for which the implementation was 



intended and we are not able to increase the H-H dis- 
tance beyond 4 A. Nevertheless, our dissociation curves 
are in good agreement with Ref. 65] and illustrates the 
dramatic failure of RPA for the atomic limit of open shell 
systems. 



IV. OUTLOOK 

In the case of adsorption of graphene on metal sur- 
faces, RPA seem to yield results that are in better agree- 
ment with experiments than both semi-local and effective 
non-local vdW density functionals. This is not surpris- 
ing since both covalent and dispersive interactions are 
important for these systems and the results seems to 
be in accordance with calculations for two-dimensional 
material a 20 ' 21 ' 35 where RPA predicts the correct inter- 
layer binding distance. However, it is well established 
that RPA does not describe covalent interactions very 
well and significantly underestimate the atomization en- 
ergies of molecules-^ and cohesive energies of solids. 18 
One would therefore expect that the dispersive interac- 
tions (far from the surface) are very well represented, 
whereas the covalent interactions (close to the surfaces) 
are less accurate. In particular, the cases of Ni(lll) and 
Co(0001) exhibit two minima which are very close in en- 
ergy (5 and 17 meV respectively). It is highly likely that 
RPA underestimate the depth of the chemisorption min- 
imum compared to the physisorption minimum and it is 
thus expected that the exact energy curves would have 
even deeper minima close to the surface. 

Since graphene on metals are being used for bench- 
marking new van der Waals functionals , 29 i 31 i 57 it is ex- 
tremely important to improve the description of such sys- 
tems beyond RPA. One line of development in this direc- 
tion is to add a second order screened exchange term to 
the RPA correlation energy^ This approach exactly can- 
cels the RPA one-electron self-correlation and improves 
molecular atomization energies slightly, but destroys the 
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accurate description of static corelation in the atomic 
limit of molecular dissociation 15 . A somewhat orthogo- 
nal line is to improve the approximation for the inter- 
acting response function within TDDFT by introducing 
an xc kerneliSr— The most sophisticated development in 
this direction is the full timc-dcpcndent EXX approach 33 
which is free of one-electron self-correlation, improves at- 
omization energies compared to RPA, and reproduces the 
correct atomic limit of static correlation. However, this 
approach may easily become prohibitly heavy due to the 
evaluation of a frequency dependent EXX kernel and it 
is not clear if the method can be generalized to periodic 
systems. On the other hand it has been shown that the 
correct dynamic properties of the xc kernel is not of vital 
importance for total energy calculations^ and one could 
simply try to use an adiabatic local xc kernel. However, 
as shown in Ref. [62j, all local kernels introduce a di- 
vergence in the pair distribution function, which makes 
it very hard to converge correlation energies and deterio- 
rates the accuracy of total energy calculations. We have 
recently shown that the divergence can be removed by a 
density dependent renormalization of adiabatic kernels, 
which defines a new class of explicit non-local adiabatic 
kernels. 34 So far this approach has been shown to sig- 
nificantly increase the accuracy of molecular atomization 
energies compared to RPA and it will very interesting to 
see if it performs equally well for periodic systems and 
graphenc on metal surfaces in particular. 
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Appendix A: Convergence of RPA calculations 

Here we will briefly discuss a few issues regarding con- 
vergence of some of the RPA calculations presented in 
this paper. 



1. Graphene on Ni(lll) 

The computational time of RPA calculations, scales as 
the number of fc-points squared, since the expression ([5]) 
involves a sum over both g-points and fc-points. Such 
scaling makes convergence with respect to fc-point sam- 
pling much more cumbersome than for standard DFT 
calculations. In particular, ab initio calculations of sys- 
tems involving graphene may require a high fc-point sam- 
pling to resolve the Dirac cone and RPA calculations of 
such systems may easily become very computationally 
demanding. Furthermore, It is not possible to perform 



an absolute convergence of the cutoff energy and extrapo- 
lation is needed in order to estimate the converged corre- 
lation energy. Since graphene on metal surfaces only bind 
by ~ 50 — 100 meV per C atom, the energy curves are 
easily destroyed by noise from the extrapolation scheme, 
which typically has an accuracy of ~ 10 meV. However, 
the extrapolation may be avoided when evaluating energy 
differences between systems of similar electronic struc- 
ture, but careful convergence tests are needed to assess 
such behavior. 

In Fig. [7] we show various convergence test for 
graphene on Ni(lll). The HF energy curves are seen 
to be highly dependent on both fc-point sampling and 
Fermi smearing. In fact, it seems extremely difficult to 
converge the fc-point sampling for the pure HF energy 
curves. Nevertheless, when the RPA energy is added the 
results are less sensitive to fc-point sampling and Fermi 
smearing and converge much more rapidly. This behavior 
is most likely due to the error cancellation between the 
expressions ([2]) and (|3]). In particular the non-interacting 
response function has been added and subtracted from 
Eq. (fl]) and should be evaluated at the same fc-point 
sampling in the two expressions. It should be remarked 
that the potential energy curve with 16 x 16 fc-point sam- 
pling is nearly identical to the one obtained in Ref. [3l[ 
with 19 x 19 fc-point sampling and we regard the en- 
ergy curve as converged. We note that the RPA energy 
curves are still more sensitive to Fermi smearing than the 
semi-local functionals and the van der Waals Functional, 
where a Fermi smearing of 0.1 eV is sufficient for a con- 
verged result. We also show the energy curves evaluated 
at cutoff energies in the range 150 — 250 eV with a 8 x 8 
fc-point sampling. The largest change is seen when in- 
creasing the cutoff from 150 eV to 200 eV. In the present 
paper we have evaluated the RPA energy curves for the 
small unit cells (Ni, Cu, and Co) using 200 eV cutoff and 
the large unit cells (Pd, Pt, Au, Ag, and Al) using 150 
eV cutoff. The choice of 150 eV for the large unit cells 
may not be quite enough for detailed convergence, but 
since the energy curves for these systems do not have 
much structure we believe that the results give the cor- 
rect qualitative features of the RPA energy curves with 
a correct equilibrium distance and binding energy which 
is within 5 meV of the converged result. 



2. Atomization energies of molecules 

In some respects, convergence of molecular atomiza- 
tion energies seems somewhat simpler than correlation 
energies of bulk systems, since one does not have to worry 
about fc-point sampling and Fermi smearing. On the 
other hand, convergence of unit cell size may become 
a major problem for a plane wave implementation. Fur- 
thermore extrapolation of cutoff energies is essential for 
molecular systems and it may be hard to obtain accurate 
results from two-point extrapolations using Eq. (|11[) . In 
general, the accuracy of the extrapolation is increased 
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FIG. 7: Energy curves for graphene on Ni(lll). Upper left: HF using 8x8 fc-point sampling with different Fermi smearings. 
Upper right: HF+RPA using 8x8 fc-point sampling and E^ A = 150 eV with different Fermi smearings. Middle left: HF 
at different fc-point samplings. Middle right: HF+RPA at different fc-point samplings using and E? u p t A = 200 eV. Bottom: 
HF+RPA at different cutoff energies for the RPA correlation energy using 8x8 fc-point sampling 



with increasing unit cell size, since the increased num- 
ber of plane waves at a given cutoff energy results in a 
smoother cutoff dependence. 

Here we show a few convergence tests exemplified by 
the atomization energy of the CO molecule. In Fig. [8] 
we show the correlation energy contribution to the atom- 
ization energy at different unit cell sizes and the extrap- 
olated results. The extrapolated results were obtained 
by a two-point extrapolation using Eq. pip with two 
subsequent cutoff points. In the left column, slightly dif- 
ferent unit cells have been used in the evaluation of O, C, 



and CO correlation energies and the atomization energy 
converges rather roughly. In the right column, the same 
unit cell was used for O, C, and CO and a much smoother 
convergence is observed. When evaluating energy differ- 
ences one can thus benefit from error cancellation when 
the same unit cell size and set of plane waves are used. 
In the present case, we see that the results are converged 
to within 20 meV when the distance between periodic 
images exceeds 6 A. However, for energy differences orig- 
inating from different unit cell sizes, the error from the 
extrapolation is larger than 0.1 eV. For the same reason, 
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it is much harder to obtain a high accuracy when eval- take energy differences between identical unit cells, 
uating the cohesive energies of solids where one cannot 
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